MISSISSIPPI STATE* UNIVERSITY / National Science Foundation 


FINAL REPORT 

NAG-1 -1061 


/AJ 



(\6f. 


Enhancement Of Surface Definition 
And Gridding In The Eagle Code 


ENGINEERING 3 
RESEARCH CENTER a 

COMPUTATIONAL 
FIELD SIMULATION 

COMPLEX GEOMETRY/ COMPLEX PHYSICS 



Dr. Joe F. Thompson, Principal Investigator 
Brian A. Jean, Graduate Research Assistant 
P.O.Box 6176 

Mississippi State, MS 39762 


{ NASA-CR- 18944 5 ) ENHANCEMENT OF SURFACE 
DEFINITION AND GRIDDING IN THE EAGLE CODE 
Final Technical Report, Aug. 1989 - Sep. 

1991 ( Mi ss i ss i pp i State Univ.) 95 p 

( CSC L 098 G3/61 


N92-14607 

Unci as 
0053187 


SECURITY CLASSIFICATION OF THIS PAGE 


Va. REPORT SECURITY CLASSIFICATION 

unclassified 


2a. SECURITY CLASSIFICATION AUTHORITY 


REPORT DOCUMENTATION PAGE 


lb. RESTRICTIVE MARKINGS 
N/A 


2b. DECLASSIFICATION / DOWNGRADING SCHEDULE 

N/A 


4. PERFORMING ORGANIZATION REPORT NUMBER(S) 


3. DISTRIBUTION /AVAILABILITY OF REPORT 


S. MONITORING ORGANIZATION REPORT NUMBER(S) 


6a. NAME OF PERFORMING ORGANIZATION 6b OFFICE SYMBOL 7a. NAME OF MONITORING ORGANIZATION 
1 Engineering Research Center (if applicable) 

Mississippi State University 


' 6c. ADDRESS {City. State, and ZIP Code) 

: Box 6176 

Mississippi State, MS 39762 


^3a. NAME OF FUNDING /SPONSORING 
“ ORGANIZATION 

j Langley Research Center, NASA 


; 8c ADDRESS {City. State, and ZIP Code) 

NASA Langley Research Center 
SSD M/S 366 

Hampton, VA 23665-5225 


7b. ADDRESS (City. State, and ZIP Code) 


8b. OFFICE SYMBOL I 9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 
(If applicable) 1 


10. SOURCE OF FUNDING NUMBERS 


PROGRAM I PROJECT 

ELEMENT NO. I NO. 



11. TITLE (Include Security Classification) 

Enhancement of Surface Definition and Griddin; 

g in the EAGLE Code 

12. PERSONAL AUTHOR(S) 

I Dr. Joe F. Thompson 

Li 3S. TYPE OF REPORT j 

[ final technical \ 

[l 3b. TIME COVERED ! 

FROM Aug 89 TO Sept 9 ll 

1 1 4. DATE Df REPORT (Year. Month. 0*v\ <15 PAGF COUNT 

! 27 , November 1991 1 94 


16. SUPPLEMENTARY NOTATION 



COSAT1 CODES 


GROUP I SUB-GROUP 



18. SUBJECT TERMS (Continue on reverse if necessary and identify by block number) 


19. ABSTRACT (Continue on reverse if necessary and identify by block number) 

Algorithms for smoothing of curves and surfaces for the EAGLE grid generation 
program are presented. The method uses- an existing automated technique which 
► detects undesirable geometric characteristics by using a local fairness criterion. 

The geometric entity is then smoothed by repeated removal and insertion of spline 
knots in the vicinity of the geometric irregularity. The smoothing algorithm is 
formulated for use with curves in B-spline form and tensor product B-spline surfaces. 


20. DISTRIBUTION /AVAILABILITY OF ABSTRACT 21. ABSTRACT SECURITY CLASSIFICATION 

□ UNCLASSIFIED/UNLIMITED □ SAME AS RPT. □ DTIC USERS 


22a. NAME OF RESPONSIBLE INDIVIDUAL 22b. TELEPHONE (Include Area Code) 22 c. OFFICE SYMBOL 

Joe F. Thompson, PhD (601)325-7299 


OD FORM 1 473, 84 MAR 83 APR edition may be used until exhausted. SECURITY CLASSIFICATION OF THIS PAGE 

All other editions are obsolete. " ~ 






















MISSISSIPPI STATE U N I V E R S I T Y / N a t i o n a I Science Foundation 


FINAL REPORT 

NAG- 1-1061 


Enhancement of Surface Definition 
and Gridding in the EAGLE Code 





TABLE OF CONTENTS 


DEDICATION 


ii 


ACKNOWLEDGMENTS . iii 

LIST OF FIGURES vi 

CHAPTER I 

INTRODUCTION 1 

CHAPTER II 

CURVES 5 

The Parametric Cubic Curve 6 

Bezier Curves 8 

The de Casteljau Algorithm 8 

Bernstein Polynomials 10 

Bezier Curves Using Bernstein Polynomials 13 

Derivatives of Bezier Curves 14 

Subdivision of a Bezier Curve 16 

Parametric Spline Curves 17 

Global and Local Parameters 18 

The Cubic Polynomial Spline 18 

Composite Bezier Curves 21 

Continuity at Junction Points 21 

Cl Continuity 22 

C2 Continuity 23 

B-spline Curves 24 


Quadratic B-spline Curves from Composite Quadratic Bezier Curves 


24 

C2 Cubic B-spline Curve from Composite Cubic Bezier Curve 25 


The B-Spline Basis 28 

Derivatives of B-spline Curves 30 


CHAPTER III 

SURFACES 33 

Bicubic Hermite Surface Patches 33 


iv 

PRECEDING PAGE BLANK NOT FILMED 


Tensor Product Bezier Surfaces 34 

Tensor Product B— spline Surfaces 36 

CHAPTER IV 

INTERROGATION AND SMOOTHING 39 

Interrogation of B-spline Curves 39 

Interrogation of Surfaces 40 

Smoothing of B— spline Curves 42 

Smoothing of Tfensor Product B— spline Surfaces 46 

CHAPTER V 

RESULTS AND CONCLUSIONS 47 

Curve Smoothing Results 47 

Waisted Body 47 

Digitized Curve 48 

NACA 0012 Airfoil 48 

Surface Smoothing Results 49 

Flat Plate 49 

Sculpted Surface 60 

F-15 Aft Quarter 60 

Conclusions and Recommendations 60 

REFERENCES 86 


PRECEDING PAGE BLANK NOT FILMED 


LIST OF FIGURES 


Figure 2.1 The de Casteljau algorithm for n=3 9 

Figure 2.2 Interplay between global and local parameters for a composite 
curve 19 

Figure 2.3 The Cl condition for composite Bezier curves 23 

Figure 2.4 The C2 condition for two degree— n Bezier curves 24 

Figure 2.5 Quadratic B-spline curve with its Bezier control polygon. ... 26 

Figure 2.6 Cubic B— spline curve with its control vertices and Bezier points 
27 

Figure 3.1 Parametric space of a bicubic Hermite surface patch 34 

Figure 3.2 Surface created by a moving and deforming curve 35 

Figure 3.3 Knot matrix of a surface spline which can be represented by a 
tensor product 37 

Figure 3.4 Knot matrix of a surface spline which cannot be represented by a 
tensor product 38 

Figure 4.1 Geometric interpretation of the B-spline smoother. 44 

Figure 5.1 B-spline control polygon for the waisted body geometry before and 
after smoothing 52 

Figure 5.2 Curvature plots before and after fairing for the waisted body . 53 

Figure 6.3 Effects of smoothing on the metric coefficient rjy for the waisted 
body 64 


Figure 5.4 Magnitude of total point displacement on the waisted body due to 
smoothing 55 

Figure 5.5 Original digitized curve data 56 

Figure 6.6 Results of smoothing with interrogation and point picking ... 67 

Figure 6.7 Curvature plots for unsmoothed data and data smoothed using the 
interrogation algorithm 68 

Figure 6.8 Results of smoothing without using interrogation and point 

picking 69 

Figure 6.9 Curvature plots for original data and data smoothed without 

interrogation and picking 60 

Figure 6.10 B-spline control polygon for original and smoothed NACA 0012 
airfoil 61 


Figure 6.11 Closeup of leading edge of the NACA 0012 airfoil showing effects 
of smoothing on the shape of the curve 62 

Figure 6.12 Curvature plots for the original data and smoothed data from the 
NACA 0012 wing section 63 

Figure 6.13 Cp vs. X/C plot for original NACA 0012 data 64 

Figure 6.14 Cp vs. X/C plot for the smoothed NACA 0012 data 66 

Figure 6.16 Pressure field near the leading edge of the original NACA 0012 
airfoil at zero degrees A.O.A. 66 

Figure 6.16 Pressure field near the leading edge of the smoothed NACA 0012 
airfoil at zero degrees A.O.A. 67 

Figure 6.17 Pressure field near the leading edge of the original NACA 0012 
airfoil at five degrees AO.A 68 

Figure 6.18 Pressure field near the leading edge of the smoothed NACA 0012 
airfoil at five degrees A.O.A 69 

Figure 6.19 Planer surface with geometric irregularities 70 


Figure 6.20 Color contours of absolute curvature on the perturbed planer sur- 
face 71 


Figure 5.21 Planer surface after smoothing 72 

Figure 6.22 Color contours of absolute curvature on the smoothed planer 
surface 73 

Figure 5.23 Original surface data for the lifting body surface patch 74 

Figure 5.24 Absolute curvature on the original lifting body surface spline 75 

Figure 6.26 Smoothed data for the lifting body surface patch 76 

Figure 6.26 Absolute curvature on the smoothed lifting body surface spline 

77 

Figure 6.27 First principle curvature on the smoothed lifting body surface 
spline 78 

Figure 6.28 Gaussian curvature on the smoothed lifting body surface spline 
79 

Figure 6.29 Original surface data from the digitized F-15 aft quarter ... 80 

Figure 6.30 Gaussian curvature on the original F-16 aft quarter surface 

spline 81 

Figure 6.31 Absolute curvature on the original F— 15 aft quarter surface 

spline 82 

Figure 6.32 Smoothed F— 16 aft quarter surface data 83 

Figure 6.33 Gaussian curvature on the smoothed F— 16 aft quarter surface 
spline 84 

Figure 5.34 Absolute curvature on the smoothed F— 16 aft quarter surface 
spline 85 


CHAPTER I 


INTRODUCTION 

Recent years have seen a dramatic increase in the use of computational 
field simulation (CFS) for the solution of physical field problems which can be 
modeled by systems of partial differential equations (PDEs). At present, com- 
putational fluid dynamics (CFD) is the most widely applied and most thor- 
oughly understood area of CFS. Current applications in CFD include flow so- 
lutions about bodies composed of complex geometrical shapes including 
sculpted surfaces. Correct modeling of physical phenomena in CFD requires 
the use of an accurate description of the geometry involved in the problem. 
Therefore, discretized curves and surfaces representing the geometry of a 
physical field problem must be as accurate as possible. 

Geometry data for CFD problems are often given as a set of digitized 
data points or as surface patches generated by computer aided drafting and 
design (CADD) software. In either case, the geometrical entities of interest 
are seldom analytical curves or surfaces. Digitized data sets usually contain 
errors resulting from the digitizing process. These errors manifest themselves 
as dimples, bumps, or ridges on the surface. CADD data bases may also con- 
tain errors due to lack of second or third derivative continuity at surface patch 
interfaces, again resulting in bumps or ridges. These errors are usually quite 
small and are not easily detected unless the surface or curve is plotted at full 
scale. However, it is not always practical or possible to examine a full scale 
plot or model since producing them is both expensive and time consuming. 


Therefore, a method of interrogation is needed which will detect small flaws 
without using a full scale representation of the geometry involved. After an 
irregularity has been detected, it is then desirable to develop a means to cor- 
rect the error without significantly perturbing the original geometry. 

In order to effectively interrogate a geometric entity, a quantity describ- 
ing the qualities sought must be defined. This quantity can then be calcu- 
lated at various locations on the geometry and its quality assessed. Several 
means of interrogation have been suggested. For detection of irregularities in 
surfaces, Kaufinann and Klass [1] use reflection lines. Beck [2] uses a variety 
of methods including contour plots, shaded images, and maps of principal cur- 
vature. Hoschek [3] uses k-orthotomic curves for interrogation of planar 
curves, and Renz [4] uses second divided differences. Farin suggests the use of 
curvature plots (Refs. [6], [6], [7]) since they are highly sensitive to changes in 
curve shape and they allow easy detection and location of inflection points. 
However, signed curvature applies only to planar curves and is, by definition, 
nonnegative for 6pace curves; as a more general approach, third derivative 
may be used (Ref [8]). 

Removal of irregularities in a curve or surface basically involves 
smoothing or fairing the geometry. Many smoothing techniques have been put 
forth. Smoothing of Bezier curves and surfaces is discussed by Hoschek[3]. 
Kjellander provides a method for smoothing cubic spline curves in Reft9], and 
bicubic parametric surface patches are discussed by Kjellander in RefflO]. A 
method for curves composed of digitized data points is given by Renz[4]. Sev- 
eral methods for smoothing B-spline curves are presented by Klass[l], and by 
Farin[6],[6],[7],[8],[ll]. It should be noted that curve and surface smoothers 
depend on the geometry model used. 


Many mathematical methods of geometry description are available and 
they all fall into two major categories, interpolation methods and approxima- 
tion methods. Interpolation techniques are characterized by the generated 
curve or surface passing through the physical data points which describe the 
object. Interpolation schemes include Lagrange and Hermite interpolation 
[7], [12], [13], [14], [15], Ferguson patches[13], Coon’s patches[13], and transfi- 
nite interpolation[14], as well as polynomial splines[7],[12],[13],[16] and ten- 
sion splines[12]. Approximation techniques are somewhat different. Instead 
of using points that lie on the object in question, they use a set of control 
points which, in general, do not lie on the curve or surface being described. 
This may lead to some misunderstanding. It should be noted that approxima- 
tion techniques represent geometry as accurately as do interpolation methods 
( in some cases, more accurately ). The difference lies in the fact that these 
schemes store the object in question using a set of control points which do not 
lie on the object. The two most widely used approximation techniques are Bez- 
ier curves and surfaces and B— splines. Excellent treatment of Bezier methods 
is given in Ref[7]; other works on Bezier methods include 
Refs[12],pL3],[15],[16]. B-splines are covered in Refs[7],[12],[17],[16]. 

Smoothing schemes based on polynomial splines or Bezier curves re- 
quire that the entire curve or surface be modified. This is due to the fact that 
polynomial splines do not exhibit the property of local control-, Bezier curves 
do not have this property either. However B-splines do exhibit local control 
and allow the development of fairing algorithms which modify a curve or sur- 
face only in the region immediately surrounding the affected control points. 
Because B-splines exhibit the property of local control, and thus allow local 
smoothing algorithms to be developed, the geometry model used here will be 



based on the B-spline representation. B— splines have already been incorpo- 
rated into the EAGLE code (RefI18]), thus making implementation of a B- 
spline smoother a logical extension of existing capabilities. 

This work is organized as follows: A brief discussion of various curve 
modeling techniques is given in Chapter II. Chapter III is devoted to surface 
generation, and Chapter IV presents the interrogation and smoothing algo- 
rithms. Chapter V contains results obtained using the interrogation and 
smoothing routines. 



CHAPTER II 


CURVES 

Many different kinds of geometry models exist. The two major catego- 
ries are interpolation methods and approximation methods. Within each of 
those categories are methods which produce single curve or surface patches, 
as well as schemes which result in a piecewise description of the geometry. 
While methods resulting in single curve or surface patches are simple and pro- 
duce smoothly varying geometric entities, they lack the flexibility necessary to 
describe complex shapes. Polynomial segments of the degree necessary for de- 
scription of these complex geometries exhibit unwanted oscillations. High de- 
gree Bezier techniques eliminate the wiggles, however they are computation- 
ally too expensive to be practical in real-world applications RefI7]. In 
addition, single Bezier curves do not have the property of local control. Piece- 
wise techniques, on the other hand, are theoretically and conceptually more 
difficult, however they allow complex shapes to be represented and are compu- 
tationally more tractable than high degree Bezier curves. Probably the most 
general of the piecewise or spline techniques is the non-uniform rational B— 
spline[7]. 

For engineering purposes, single curve segments higher than degree 
three need not be considered. This is because cubic curve segments are true 
three dimensional curves. Since engineering efforts are restricted to geome- 
tries of three dimensions, the use of higher degree curves would simply intro- 
duce unneeded complexity and computational requirements. 



The Parametric Cubic Curve 


Mortenson[16] defines a parametric cubic ( PC ) curve as a continuous, 
single-valued mathematical function of the form: 

x = x(t) y = y(t) z = z(t ) . 

The parameter t ranges from t=0 at the start of the curve to t=l at the end. In 
vector notation, the curve is denoted as P(t) where the vector P has the compo- 
nents x, y, and z in three dimensions. Derivatives of P are taken with respect 
to the parametric variable t as follows. 

P\t) = P“(t) = 

Higher derivatives are calculated in a similar fashion, however for PC curves, 
derivatives higher than third order are meaningless because they are zero. 
Algebraically, a PC curve can be expressed as a cubic polynomial in the vari- 
able t. Equation (2.1) is the algebraic form of a PC curve. Note that the coeffi- 
cients of the t terms are vectors. 

P(f) = "s* 3 + a 2* 2 + a \ f + a o (2.1) 

A more physical definition of the PC curve is offered by the geometric form. 
This form is derived by taking the first derivative of Equation (2.1) with re- 
spect to the t and solving for the tangent vectors of the curve at the ends, P*(0) 
and P^l). Then solving Equation (2.1) for the end points of the curve, P(0) 
and P(l). The results are given as Equations (2.2). If Equations (2.2) are 
solved for the algebraic coefficients a, Equations (2.3) result. Substituting 
these back into Equation (2.1) and grouping the result according to vectors P 
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and P l , results in Equation (2.4:) 

P( 0) = « 0 

P(l) = a Q + a 1 + a 2 + a 3 

( 2 . 2 ) 

P‘( 0) = fl, 

P'(l) = aj + 2a 2 + 3 a 3 


a 0 = P(0) 

«1 = no) (2.3) 

a 2 = — 3P(0) + 3P(1) - 2P'(0) - P'(l) 
a 3 = 2P(0) - 2P(1) + P'( 0) + P'(l) 

P(/) = (2Z 3 - 3 1 2 + 1)P(0) + (- 2^ + 3r 2 )P(l) + (? 3 - 2z 2 + r)P'(O) 

+ (r 3 - r 2 )P'(l) (2 .4) 

The coefficients of the P and P l terms can be interpreted as blending func- 
tions. These functions relate the physical curve to the parametric space of the 
curve. If the coefficients of P(0), P(l), P^O), and P^l) are denoted Fi, F2, F3, 
and F4 respectively, then Equation (2.4) becomes 

P(t) = FjP(O) + P 2 P( 1) + F S P‘(P) + P 4 P'( 1), (2.6) 


where 


In matrix form, 


Pit) = [f 3 


t 2 t 1 


2 

3 

0 

1 


- 2 1 

3-2 
0 1 

0 0 


1 

- 1 
0 
0 



( 2 . 6 ) 


Bezier Curves 


The de Casteljau Algorithm 

The de Casteljau algorithm generates an n— th degree polynomial curve 
by repeated linear interpolation. The algorithm is stated as follows. 


Given points bo,...,b n in 3 space and a real number t, 
b r St) = (1 - t)b\-\t) + tb r -\{i) 


r = 1 

i = 0 - r 


(2.7) 


where 


b^it) = b t . 


The final step in the algorithm ( with r=n and i=0 ) will yield the point on the 
curve at the parameter value t. Curves produced in this way are Bezier curves 
and have several important properties, some of which are discussed in Ref [7]. 
These properties include: 

» Affine invariance. 

» Invariance under affine parameter transformations. 

» Convex hull property. 

» Endpoint interpolation. 
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A graphical interpretation of the de Casteljau algorithm for the cubic 





Figure 2.1 The de Casteljau algorithm for n=3. 

Listed below are two examples of expansions of the de Casteljau algorithm for 
the quadratic and the cubic cases. 

n = 2 

b l 0 (t ) = (1 - t)b Q + tb x 
Z»](r) = (1 - t)b 1 + tb 2 

b 2 0 (t) = (1 - t)b l Q + tb\ 

Substituting the first two Equations into the second, a closed formula results: 
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b 2 0 {t) = (1 


— t) 2 b 0 + 2r(l — t)b j + t 2 b 2 


n - 3 

* 1(0 = (1 ~ 0 * o + tb \ 
b\{t) = (1 — r)£>j + 

p 2 it) ~ ( 1 . — 0*2 
*1(0 = (1 - 0*1 + tb\ 

*?(0 = (1 - 0 *] + f*J 
*1(0 = (1 - 0*1 + **? ' 

Performing a substitution process similar to that for the quadratic case results 
in Equation (2.8): 

blit) = (1 - 0 3 * o + 3f(l - 0 2 *i + 3r 2 (l - 0*2 + ' 3 * 3 (2.8) 


Bernstein Polynomials 


Bernstein polynomials are defined explicitly by 


5”(0 = C(n,00(l - 0" _i ' 


Where 


t e [0 . l] 


C(n,i) = 



n\ 


J i\(n — /)! 


0 


, z € [ 0,/z ] 

, otherwise. 


(2.9) 




A recursion formula for the Bernstein polynomials can be derived as follows 
and the result is given as Equation (2.10): 


B i (0 = (?)*1 - t r~ i 

= (” 7 1 ) r ( i ~ / ) n_I + ~ i ) r *(i - r )" _/ 

= ( i - 1 )B n i -\ty + 1 te to. i] (2.io) 

The results of evaluating Equation (2.9) for the cases n=2 and n=3 are 
given by Equations (2.11) and (2.12) respectively. These correspond to the 
quadratic case and cubic case respectively: 

= ( 1 - n 2 

B\(t) = 2/(1—/) 

Bjit) = f 

B 2 0 (t ) = ( 1 - / ) 3 
Blit) = 3/ ( 1 - / ) 2 
B 2 0 it) = 3/2(1-/) 

B 2 Q {t) = Z 3 

Some properties of Bernstein polynomials, from Refs [19], [20], [21], include 
the partition of unity, the positive property, endpoint interpolation, the recur- 
sive property, and symmetry: 

» Partition of Unity: see Equation (2.13). 

» Positive Property: For 0 < t < 1 Bernstein polynomials are positive. 

» Endpoint Interpolation 
» Recursive Property: see Equation (2.10). 


( 2 . 11 ) 


( 2 . 12 ) 


5-1*. 


(2.13) 


I>"« = 1 

£ = 0 

Proof: 

n n 

i = a+d- w = e 7Vo - = E w 

i=0 ' ■ 1 = 0 


fi”(G) = 1 iff i = 0 
B n t { 1) =1 iff i = « 

Symmetry: From the symmetry property of the binomial coefficients, 


C( n, n- i) = C{ n, i ), 

t ) = C( n, n — i )*»"'( 1 - f V 

= C( »,/ )( 1 - 0‘t 1 - ( 1 - O l" -1 ' 

= 1 “ 1 ) • (2.14) 

Derivatives of Bernstein polynomials are determined as follows: 


| fi"(o = i an ,,: ) /»(!-/ )“-■ 


= C(n,i)[ < 1 - 1 - ( n - / ) !"( 1 - / ) 


n~i— 1 


] 


i n\ 


i! (n- i )! 

n ( n — 1 )! 
(>-!)!(«-/)! 


[<'-'( 1 - < ] - ( ,-i ( - t ( 1 - ' r - '' 1 ] 

( <'-( i - ' « , n-' r-'- ) 


5"(0 = n[^r/ ( r ) ] 


AZ“ 1 


(2.15) 


Bezier Curves Using Bernstein Polynomials 
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The Bezier curve is a simple approximation technique which takes the 
form of Equation (2.16): 

i— n 

P{t) = ^ bpy) t E [0, 1] (2.16) 

/ = 0 


Here the bj are the control vertices of the curve, and the B(t) are the Bernstein 
polynomials. A word about the subscript notation in the above Equations is 
in order. Since the Bernstein polynomials depend on the number of control 
vertices used in the curve description, the double subscript notation is neces- 
sary. The subscript n denotes the degree of the curve ( i.e. n=2 for quadratic, 
n=3 for cubic, etc. ). The subscript i indicates the control vertex associated 
with the particular Bernstein polynomial. If the Bezier formulation of Equa- 
tion (2.16) is expanded for the cubic case, Equation (2.17) results: 


P(t) = (1 - t) 3 b 0 + 3r(l - t) 2 b 1 + 3^(1 - t)b 2 + t 3 b 3 


(2.17) 


Equation (2.17) was obtained directly using the Bernstein polynomials, how- 
ever, Bezier curves were first introduced in the form of a recursion formula 
( the de Casteljau algorithm ). If one compares Equation (2.17) with Equation 
(2.8), which was obtained from the de Casteljau algorithm with n=3, they are 
the same. If Equation (2.17) is cast in the form of a PC curve, the blending 
functions F now become 


Fi = (1 - tf 
F 2 = 3r(l - t ) 2 
F s = 3^(1 - t) 
F 4 = t 3 


and Equation (2.17) becomes 


Pit) = Fj6 0 + F 2 b 1 + F 3 b 2 + F 4 b 3 , 
or in matrix form, 

Pit) = [t 3 t 2 t 1 


1 

3 

3 

1 


3 

6 

3 

0 


3 1 
3 0 
0 0 
0 0 


- 

1 o 
<5 


*1 


b 2 

- 

cn 1 

i 


(2.18) 


(2.19) 


Because the Bernstein polynomials form a basis for Bezier curves, Bezi- 
er curves inherit several properties intrinsic to the Bernstein polynomials. 
The following is a partial list. Farin[7] provides a more complete listing. 

» Affine invariance. Because of Equation (2.13), Bezier curves are in- 
variant under affine maps. 

» Invariance under affine parameter transformations. 

» Convex hull property. Because of the positive property of the Bern- 
stein polynomials. 

» Endpoint interpolation. 

» Symmetry. Due to Equation (2.14). Algebraically, this is stated as 
Equation (2.20): 




i= 0 


i ‘=0 


f»?( 1 - t) 


( 2 . 20 ) 


Derivatives of Bezier Curves 

From Equations (2.15) and (2.16) the derivative of a Bezier curve can be 


determined: 


15 




j t b'\ /) = »£( B i-i^ o - B i~ l ( o ) b i 

i=0 

This can be simplified to yield 

rt — 1 
1 = 0 

Where 

Ai >i = *,+ i - 


( 2 . 21 ) 


Repeated application of Equation (2.21) produces a general formula for deter- 
mining higher order derivatives of a Bezier curve. First, the forward differ- 
ence operator 6 is generalized in Equation (2.22). The formula for the r-th 
derivative of a Bezier curve is given in Equation (2.23). 


( 2 . 22 ) 

(2.23) 

Of particular interest are the derivatives at the ends of the curve. These cor- 
respond to t=0 and t=l. Evaluating Equation (2.23) at t=0 and t=l, and recal- 
ling the endpoint interpolation property of Bernstein polynomials, Equations 
(2.24) and (2.25) result: 


A r b,= J C(r,k)( - 1 V~ k b i+k 


k= 1 


ri i 

b n (t) = 7 ” ! Y A r b, B r }~ r ( t ) 

dt r (n - r)\ A-, J j v ' 


7'=o 
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d r 
dt r 


d r 
dt r 


b\ 0 ) = 
b n ( 1 ) = 


>V. , ru 

/ \,Zl 

( n — r )! 0 

/j! ^ ri 

^ bn — r 


( «-r)! 


(2.24) 

(2.26) 


Subdivision of a Bezier Curve 

Subdivision is the process by which a single Bezier curve of degree n, 
defined over the interval t£ [0,1] is divided into two curves of the same degree 
defined over the intervals t£ [0,a] and t£ [a,l]. The following development is 
due to Farin [7]. 

The first curve segment (t=0tot=a) will be denoted as a 11 and will 
have control vertices ao,...,a n . A local parameter v will be introduced such that 
v = t/a; therefore, as v varies from 0 to 1, t varies from 0 to a. Since the curve 
segment a 11 is part of the same curve segment as b 11 , all of their derivatives at t 
= v = 0 must agree: 


jL„n (0) = AL a «( 0); r = 0,...,n 


dv 


(2.26) 


From Equation (2.23), the derivatives in Equation (2.26) depend only on the 
control points ao,...,a r and bo,...,b r Now, if these two sets of r+1 points are tak- 
en as control points of two degree— r Bezier curves, then then two curves are 
the same for all v and t because of Equation (2.26): 


«o( v ) = b r 0 (t ) 


Finally, from Equation (2.27) with v=l and t=a, 


(2.27) 


The left hand side of Equation (2.28) is actually a r ( the control points of the 
first Bezier curve segment ). Therefore, 

aj = V Q {v) . (2.29) 

Equation (2.29) is known as the subdivision formula for Bezier curves. This 
allows the de Casteljau algorithm to compute the control polygon for the curve 
segment corresponding to the interval [0,a]. Because of the symmetry property 
given by Equation (2.20), the polygon of the curve segment corresponding to 

the interval [a, 1] is composed of the vertices . 

Parametric Spline Curves 

The spline curve derives its name from a drafting tool called a spline. 
This is a thin piece of wood, metal, or some other flexible material which is 
made to pass through a series of established control points on the drawing of 
interest. While the spline is held in place by weights, the designer uses it as a 
guide to draw a smooth curve passing through the designated points. 

The mathematical spline is a smooth piecewise polynomial curve 
which can be drawn through any set of points. The degree of the spline is de- 
termined by the degree of the curves which make up each span of the spline. 
There should be no kinks or rapid changes in curvature in the spline curve. 
The spline curve will be denoted by s(u) where u is the global parameter of the 
spline. On each span of the spline, the curve will be written in terms a local 


parameter t. Therefore, on an particular span, the spline will be a function of 
the local parameter t so that s(u) can be written as 


x = x(t) y = y(t) z = z(r) . 

The locations in parametric space where the spans of the spline join are 
known as knots. Derivative continuity is prescribed at the knots. Quadratic 
splines are only capable of first derivative continuity ( C 1 ), while cubic splines 
may be C 2 . The set of knot values for the entire spline is called the knot vector 
and will be denoted by U. Individual components of the knot vector are given 
by Uj where i denotes the location of the knot within the vector. 


Global and Local Parameters 

It will be useful to define two parameters for a spline curve, a global 
parameter u and a local parameter t. The global parameter u varies continu- 
ously over the entire spline curve, while the local parameter t is constrained to 
the interval [0,1] on each span of the spline. The relation between the global 
and local parameters is illustrated by Figure 2.2 and Equation (2.30). 


t = 


u — u i 
u i + 1 " u i 


(2.30) 


It is easy to confirm that over any parametric interval [uj,ui + i], t will vary over 
the interval [0,1]. Thus, the spline curve can now be written as Sj(t) where the 
i denotes the i-th segment of the spline, and s(u) = Sj(t). 
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Ui_i uj Ui +1 ui+2 


Figure 2.2 Interplay between global and local parameters for a composite 
curve 


The Cubic Polynomial Spline 

The cubic polynomial spline splines a set of data points, xj, with the pa- 
rameter u. On any interval [uj,ui+i], the spline curve s(u) = Si(t) is given by 


sfr) = k<i - ty 


+ l c i + 1 


t 3 + ( Xi - ic/)(! - 1) + (*/+i - k+ 1)* • 


( 2 . 31 ) 
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The Cj in Equation (2.31) are found by the tridiagonal system given below in 
Equation (2.32). 

c i _ 1 + 4c i + c i+ ! = 6(x i+ j - 2x f + ,) . (2.32) 

Since the solution of Equation (2.32) is necessary to compute the spline, any 
change in the data points X{ results in changing the entire spline curve. This 
means that the polynomial spline lacks the properly of local control. 

The ends of the spline require some special treatment. Several types of 
end conditions can be prescribed including zero curvature or natural end con- 
ditions, constant curvature or quadratic end conditions, and specified slopes. 
Details on the various types of end conditions can be found in Reft 22], 

Perhaps the most widely used end condition type is the quadratic end 
condition. It will be used here as an example. In this case, Equation (2.32) is 
replaced by 

c 2 ~ x 3~ 2 x 2 + X 1 

C NI— 1 = X NI ~ ^ X NI- 1 + x NI-2 

for i = 2 and i = NI where NI is the number of data points on the spline curve. 
The values of c for the first and last points are determined by extrapolation as 
follows: 

Cj ' — 2c2 C3 

c Nl ~ - C N! _ 2 


Composite Bezier Curves 


As mentioned earlier, single Bezier curves are not capable of describing 
complex shapes unless they are of prohibitively high degree. However, the use 
of composite Bezier curves will allow complex shapes to be modeled while 
keeping computational requirements within reason. In addition, the concepts 
of quadratic and cubic B-spline curves can be developed directly from compos- 
ite Bezier curves. The following development is based on that of Farm [7] and 
uses his notation. 


Continuity at Junction Points 

Since individual Bezier curves are polynomials, they are continuous in 
all derivatives. However, for composite Bezier curves, continuity at the junc- 
tion points must be prescribed. The cases of C 1 and C 2 continuity conditions 
at these junction points provide a basis for the development of quadratic and 
cubic B-spline curves respectively. 

Let two Bezier curves s 0 and Si be described by control polygons 
bo, —.bn and b n ,...,b2 n respectively and defined over the intervals [uq.ui] and 
[ui,U 2 h If these curves are the product of a subdivision process, this means 
that the curves s 0 and Si are part of a single polynomial curve defined over the 
interval [uo,U 2 ]. Therefore the control vertices of the two curves are related by 
Equation (2.33): 


b n +i = i = 0,...,n 


(2.33) 


where 
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Since the two n— th degree Bezier curves are part of one n— th degree 
polynomial, the two curves must agree in all derivatives up to the n-th deriva- 
tive. From Equation (2.23), the n— th derivative of a Bezier curve depends only 
on the surrounding n control vertices. Therefore, using that knowledge and 
Equation (2.33), the C r condition for the junction point of two Bezier curves 
can be constructed. 


Given two Bezier curves defined over the intervals [u 0 ,ui] and 
[ui,U 2 l and by the control polygons bo,...,b n and b n ,...,b2 n respec- 
tively. They are r times continuously differentiable at ui iff 


b n+i = i = 0,...,r 


where 


t 


u. 



o 

0 


(2.34) 


The C r condition can also be formulated by equating derivatives at the junc- 
tion point. From Equation (2.23), 

Ai b n -i = [j~j Albnt 1 = °’ -’ r ' (2 ‘ 36) 



Note that the derivatives here are with respect to the global parameter u and 
not the local parameter t, therefore the chain rule applies. 


C 1 Continuity 

® From Equation (2.35) with i = 0,1, 

bn = b, i 

and 

• A l Ab n - l = A^b n (2.36) 

or, expanding the second equation, 


(u 2 - u x )(b n - b n _ j) = - u 0 )(b. n+ j - b n ) . 


As can be seen from Equation (2.36), the three points b n _i, b n , b n+ i must be 
colinear; however, this is not sufficient to guarantee C 1 continuity. The addi- 
tional requirement is that the three points be in the ratio &o : &i- This is illus- 
trated in Figure 2.3 . The C 1 condition can be used to.determine a formula for 


^n-l AO : ai b n+1 



Figure 2.3 The C 1 condition for composite Bezier 
curves. 

computing the location of the junction point bn given the two surrounding con- 
trol vertices b n _i and b n+ i. This formula is given by Equation (2.37) and will 
be used shortly to develop the concept of a quadratic B-spline curve. 


b n = 


'0 


fli b + 

+ A l bn ~ 1+ A 


'o 


+ A 1 


'n+l 


(2.37) 


or 


b 


n 


u 2 U \ u , U \ U 0 L 

U 2 ~ U 0° n ~ i ^ U 2 ~ u o° n+1 


C 2 Continuity 

If a given curve is C 1 at the knot ui then from Equation (2.35), the addi- 
tional requirement for C 2 continuity is that an auxiliary point d be uniquely 
defined by Equations (2.38) and (2.39). 
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b n - 1 — (1 ? i)^n-2 

^n+1 = (1 — + t\b n+ 2 

where 

_ ~ 

1 1 U 2 — Uq- 


(2.38) 

(2.39) 


Note that ti is actually the local parameter of ui with respect to the interval 
[110,112]. The polygon b n _2, d, b n+ 2 describes a single quadratic polynomial over 
the interval [uo,U2l. Figure 2.4 shows a curve that is C 2 at the junction point 
b n : 




Quadratic B-spline Curves from Composite Quadratic Bezier Curves 
If some quadratic spline curve s defined over the interval [uo,ujJ is as- 
sumed to be C 1 then the spline curve can be completely defined by the knot 
vector u and the inner Bezier points bo, bj, b3 b2i+i,...,b2L_i, b2L- The addi- 



tional points (junction points ) required to define the composite quadratic Bez- 
ier curve can be determined by the C 1 condition given in Equation (2.36). 
Solving Equation (2.36) for the junction point b 2 i in terms of the inner Bezier 
points, Equation (2.40) results: 


b n A 


£? i ft -j- ^i-\ Ij 


(2.40) 


The polygon bo, bi, b 3 ,..., b 2 i+i b 2 L_i, b2L can now be called the B- 

spline control polygon of the spline curve s. This polygon, along the its knot 
vector u, completely describe the quadratic spline curve s. To indicate that the 
control polygon now denotes a B-spline curve, the control polygon will be de- 
noted by d_i, do, di,...,dj v _i, di,. Figure 2.6 illustrates the relationship be- 
tween the B-spline control polygon and the Bezier control net for the quadrat- 
ic case. Note that the polygon for the curve is stored as follows: 



With the inner vertices given by 
d[ — b 2i + 1 i = 0, ... L — 1 . 


C 2 Cubic B-spline Curve from Composite Cubic Bezier Curve 
Consider a cubic spline curve defined over the knot sequence [uo,ul]. In 
order for the curve to be C 1 , Equation (2.40) must hold. For the cubic case, the 
C 1 condition is given by Equation (2.41). 
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The additional condition for C 2 continuity, previously given as Equations 
(2.38) and (2.39), is now given by Equations (2.42) and (2.43): 



The ends of the spline are treated such that the first and last control vertices 
of the spline are coincident with the first and last points of the curve respec- 


tively. Equations (2.41) through (2.43) are valid for i=l,...,L— 1. For the end 
points of the curve, Equations (2.44) are used. 


b 0 d_ j 


b i d Q 

b = d 0 + -r^V 




l. _ "L-l j , a L - 2 

» 3 L -2 “ iJ L _ 2 + J + A + J 


(2.44) 


'L- 1 


1 L — 2 


L-l 


^3L-1 — d L 


&3 L ~ ^i+1 



Figure 2.6 Cubic B-spline curve with its control vertices and Bezier 
points. 



The B— Spline Basis 


Further theoretical development of B— splines requires defining them in 
more analytical terms. Equation (2.45) is the Cox, de Boor recursion for B- 
spline curves. This formula can be used to calculate B— spline basis functions 
for splines of any degree. When Equation (2.45) is expanded for the cubic 


u - u i _ l 


A J n (u\ = I Z Lz i N n ~ 1 (u) + — N n ~ l (u ) 

U i + n- 1 - U i - 1 * W U i+n ~ *+» W 


U i + n ~ U 


7 / 1 - 1 / 


with 

1 if u i _ 1 < u < u L 
0 else 



(2.45) 


case. Equation (2.46) results. Using Equation (2.45), the equation of a B~ 
spline curve for a given set of control vertices and a given knot vector is given 
as Equation (2.47) where L is the number of interval on the spline. Note that 
Equations (2.45) and (2.47) are given in terms of the global parameter u of the 
spline. If these equations are cast in terms of the local coordinate t of the 
curve instead of the global coordinate u, then Equation (2.47) becomes Equa- 
tion (2.48) where t is the local coordinate of the ith segment of the curve and is 
given by Equation (2.30). For the case of a cubic curve, the matrix form of 
Equation (2.48) is given by Equation (2.49). Where M 3 is a 3x3 matrix given 
by Equation (2.60). 
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0 “ - U i ) ? r r , 

(«i+3 - “/) ("f+2 - "i) («i+l - «i) 


(« - 


( w ,+ 3 “ w i) K + 2 ~ "/) ( M / + l “ “/) 


« e [w i+1 ,« l+2 ] 


+ 


(" - W J + 1 ) 3 (l/ f + 4 - M,) 


( M i + 1 “ "/) ( M , + l ~ U i+ 2> ( W / + l “ «, + ?) ("/ + ! ~ w / + 4 ) 


N*( u ) = 


(«/ + 3 “ w ) ? (", + 4 - «/) 


("/ + 3 ~ U i ) ( w / + 3 " «/) («/ + 3 - «/+2) ( M j + 3 “ w , + 4) 
+ 


( M i+4 ~ w / + l) ( W / + 4 “ w / + 2> ( M j + 4 ~ Wf + 3) 


(«i + 4 " ") 3 


("i + 4 “ W / + l) ( W ; + 4 - w , + 2> ("i+4 " w /+3) 


w e [w . .,, 1/ , +4 ] 


0 


u £ [u it U i+4 ] 


(2.46) 
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L + n— 1 


P(u)= ]T dp?(u). 


(2.47) 


i= 0 


i+n- 1 


Pi(t) = X W « • 


(2.48) 


j=i-i 


P L {t) = [l f r 2 / 3 ] M 3 


1 


*/+ 1 
h + 2 


(2.49) 


Derivatives of B— spline Curves 

For purposes of interrogation and smoothing, the derivatives of a curve 
should be calculated with respect to the curves global parameter. However, 
Equation (2.48) and, in particular, Equation (2.49) are given in terms of a local 
parameter. This requires the use of the chain rule in order to determine deriv- 
atives with respect to the global parameter u of the curve. Consider a spline 
curve Sj(t) where t is the local parameter of the curve in the interval [upui+i]. 
The first derivative of the curve with respect to the global parameter u is given 
by Equation (2.61). Higher order derivatives are given by Equations (2.62) 
and (2.63). 

For the case of cubic B-spline curves defined by Equation (2.49), first, 
second, and third derivatives are given by Equations (2.64), (2.66), and (2.66) 
respectively. 


■mu 

m n 

m 13 

™14 

m 21 

m 22 

m 23 

m 24 

m 31 

m 32 

m 33 

m 34 

™41 

™42 

m 43 

m 44 


(“« “ 

M i+l) 2 



11 ( w /+l - u i- 2) (",+1 - 1) 

m 2 1 = “ 3m 11 

m 31 = 3m lx 


13 (Mf-j- h, + 2 ) — " i+1 ) 

w = 3 (“/ + 1 ~ “i) (“i ~ “ 1 - 1 ) 

23 (W,-! - U i + 1 ) («,-_! - « ( + 1 ) 

m = 3 ^,- + i ~ w /) 2 

33 0,-j - « ;+2 ) K-i - w i+ 1 ) 

W 43 = “( W i + l“ W /) 2 


(2.60) 


" J 12 

= 

1 - ( m n + m 13 ) 

w 22 

= 

- ( ™21 + m 23 ) 

m 32 

= 

- ( m 31 + m M ) 

m 42 

= 

- ( m 41 + m 43 + m^ 

m l4 

= 

0 

m 24 

= 

0 

™34 

= 

0 



K +] - «,) 2 

m 44 


(“i ~ u i+ 3 ) (“/ " m / +2 ) 


(“i-i 




1) ( M ;-i 


+ 


1 


4 i + 2, 


l) (“. “ «< + 2 ) (“. “ «.+?) («i+2 “ «i-l) («J+2 ~ «i) 


<fe(M) = fjg#) dt = J_^C0 
Jw dt du A ■ dt 

( Ps { u ) _ 1 ^,-(0 

(Zl ,.) 2 rfr 2 

<As(m) _ 1 ^,(0 

^u 3 (zl ( .) 3 dt 3 


(2.61) 

(2.62) 


(2.63) 
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l i - 1 




Pi(t ) = ^[0 1 2r 3r 2 


M, 


1 4 


*/+ 1 
^' + 2 I 




- / I \2 


(j-) 2 [0 0 2 6r ] M 3 


‘i- i 


* 1+1 
^i+2 I 


ft- 1 


ft+r 


\d : 


/ + 2 


(2.64) 


( 2 . 66 ) 


£pP(.t) = (4) 3 [° 0 0 6]M 3 


(2.56) 


• . CHAPTER III 

SURFACES 

Curves are defined in terms of a one— dimensional parameter space with 
global parameter u and local parameter t. Surfaces, however, are defined in 
terms of a two-dimensional parameter space with global parameters u and v 
and local parameters s and t. There are many ways to represent surface 
patches, however, this work will concentrate only on the tensor product form. 

Bicubic Hermite Surface Patches 

• Surface patches of the kind found in the early versions of the EAGLE 
surface code are bicubic Hermite patches. These are given by equation (3.1): 


P(s,t) = B{s)HB r (t ) 


(3.1) 


where 


H = 


'#<0,0) r(l,0) r u { 0,0) r u (l,0)' 
K0,1) r(l,l) r u (0, 1) r u (l, 1) 
r v (0,0) r v (l,0) r„ v (0, 0) r wv ( 1 , 0) 
r v (0, 1) r v (l, 1) r uv (0, 1) r uv (l,l) 


B(a) = 


1 - 3a 2 + 2a 3 ' 
3 a 2 - 2 a 3 
a — 2 a 1 + a 3 
— a 2 + a 3 


, a = s,t 
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Here, s and t are local coordinates of the spline and the global coordinates are 
taken as the integer indices of the grid point. The parametric space of the 
spline is shown in figure Figure 3.1 . 
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Figure 3.1 Parametric space of a bicubic Hermite surface patch 
Tensor Product Bezier Surfaces 

Recent versions of the EAGLE code include the capability of generating 
tensor product Bezier surfaces (RefI183). This type of surface patch is a simple 
extension of the Bezier curve technique. The mathematical form of a tensor 
product Bezier surface patch is given by equation (3.2). 


m n 

b mjx {s, t) = 

i=0 j=0 


(3.2) 


A set of (n+l)x(m+l) control points define the surface. The surface may be 
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generated by treating each row of the control net as a degree— n Bezier curve, 
then treating each column of the result as a degree-m Bezier curve. This 
amounts to generating a surface by using a curve moving through space and 
changing its shape. The matrix form of equation (3.2)_for the cubic case is giv- 
en by equation (3.3): 


P(s,t) = SM b2 BM T bz T T 



O 

O 

^01 

*02 

*03 

B = 

*10 

*11 

*12 

*13 

^20 

*21 

*22 

*23 


^30 

*31 

*32 

*33 

S = [ 1 s 

S 2 

S 3 ] 




t t 2 / 3 ] 


0 

0 

3 

3 


O' 

0 

0 

1 


(3.3) 



Figure 3.2 Surface created by a moving and deforming curve. 
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'Itensor Product B— spline Surfaces 

Recent versions of the EAGLE code also have the capability of generat- 
ing tensor product B-spline surfaces (Ref[18]). The relationship between B— 
spline surfaces and B— spline curves is the same as that between Bezier sur- 
faces and Bezier curves. Equation (3.4) gives the mathematical form for the 
surface patch: 


m n 

(3-4) 

i=0j=0 

Note that the basis functions N are the tensor product of the set of univariate 
basis functions given by equation (2.45). The matrix form for the cubic case is 
given by equation (3.6): 

P(s,t) = sm 3 dm\t t (36) 

where D is a 4 x 4 matrix of control points and M 3 is the same as in equation 
(2.60). The vectors S and T are composed of the local parameters s and t of the 
point to be generated. 

B-spline surfaces inherit all the properties possessed by B-spline 
curves. This is due to the fact that the tensor product method can be broken 
into a series of univariate operations each of which can be interpreted as gen- 
erating a B-spline curve. Note that there are limi tations to the tensor product 
method. Perhaps the most restrictive of these is the necessity that only one 
knot vector can be prescribed for each direction on the surface. This is to say 
that all u isoparametric curves must have the same knot vector and likewise, 
all curves of constant v must be defined by a single knot sequence. Figure 3.3 
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illustrates the parametric space of a surface which is definable by a tensor 
product method, while Figure 3.4 shows the parameter space of one which 
cannot be represented by a tensor product: 



Figure 3.3 Knot matrix of a surface spline which can be represented 
by a tensor product 
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Figure 3.4 Knot matrix of a surface spline which cannot be represented 
by a tensor product 



CHAPTER IV 


INTERROGATION AND SMOOTHING 

Chapters II and III have described some curve and surface modeling 
techniques. However, if given a particular geometry, how can one access the 
quality of that geometry, in terms of derivative continuity, and if necessary, 
improve upon it? Attempts to provide a practical solution to this problem have 
resulted in a variety of techniques ( see Chapter I ). This work, however, will 
focus on techniques provided by Farin[6],[7]. 

Interrogation of B— spline Curves 

Several means of curve interrogation have been suggested. Curvature 
plots provide a great deal of information about the shape of a curve. However, 
their use is limited to planar or two dimensional geometries because curvature 
is always nonnegative for space curves. The same two-dimensional limitation 
applies to k— orthotomic curves. What is needed is a quantity which is highly 
sensitive to small perturbations in the geometry but is not limited by spatial 
dimension. The third derivative provides such a quantity and will be used for 
interrogation of space curves. For planar curves however, signed curvature 
will be used. 

Calculation of derivatives is effected by use of equations (2.54) through 
(2.66). For planar curves, curvature is calculated via equation (4.1). For de- 
termining derivative continuity, derivatives or curvature must only be calcu- 
lated at the spline knots. All points on the interior of each span of the spline 
are defined by polynomials and thus are continuous in all derivatives. 



Farin[7] suggests that a fair curve should be composed of segments with 
smoothly varying curvature ( i.e. there should be no slope discontinuities in the curva- 
ture plot of a curve ). For planar curves, this definition is sufficient, however, for 
space curves, the second derivative will be used in lieu of signed curvature. Slope 
discontinuities in second derivative are actually third derivative discontinuities. 


_ x"(u)yXu) - y"(u)x'(u) 

1(X '(U))* + (y'(u))2]3/2 


(4.1) 


The method used to locate discontinuities is to calculate the left and right hand limits 
of the third derivative at each spline knot and then compare them. This will be the 
basis for a quantity known as local fairness and leads to the following definition: 


Definition: Let x(t) be a C 2 parametric cubic piecewise curve with u as the global 
parameter. Then the local fairness e is defined as 


£ =||jf"'(lf + ) - *”'(«-) |. 


(4.2) 


Note that e is a local quantity since it may vary with the value of the parameter u. It 
is reasonable to say that the point most in need of smoothing is the point with the 
largest value of e. 

A separate routine for interrogation of curves was not produced. The 
method of interrogation presented here was incorporated into the smoothing 
routine outlined below. 


Interrogation of Surfaces 

The maximal and minimal curvatures of a surface at a point P are 
known as the principle curvatures and are denoted by ?q and %2 • These eurva- 


tures are calculated by finding the roots of equation (4.3) where the super- 
scripts on P denote differentiation with respect to the parameter indicated. 


(EG - F 2 )x 2 - (EN + GL- 2FM)x + (LN - M 2 ) = Q 
where 

E = P u . P u L = P™ ■ n 

F = pu. P v M = P™ • n 

G = P* P v N = P™ ■ n 

and 

_ _ P u x P v 

n \P“ X P v \ ' 


(4.3) 


The two major types of surface curvature are known as mean curvature H and 
Gaussian curvature K. These are given by equations (4.4) and (4.6) respec- 
tively. 


K = XjX 2 


LN- M 2 
EG-F 2 


(4.4) 


H = + *2) = 


EN + GL - 2 FM 
2 (EG-F 2 ) 


(4.6) 


While these quantities provide a great deal of information, they have limita- 
tions. Gaussian curvature is always zero for cylinders of the form given in 
equation (4.6) while mean curvature is always zero for minimal surfaces. 


c(u, v) = (1 - u)X(v) + u[X(v) + V] 


(4.6) 


Absolute curvature, however, will always detect curvature in a surface. Abso- 
lute curvature is given by equation (4.7): 

^ abs ~ hljl + IX 2 I _ (4 7) 

A subroutine for evaluation of principal, Gaussian, mean, and absolute 
curvatures was developed for the EAGLE code. The routine is capable of eva- 
luating surface qualities of either B-spline surfaces or surfaces produced by 
the original EAGLE surface spline routine, SPLSUR. In the latter case, deriv- 
atives on the surface are calculated using the SPLINT subroutine ( see Ref 
[23] ). The result is two data files, one containing points on the surface and 
the other containing values of the five curvatures at those points. These files 
are written in formatted PLOT3D single grid form. 

Smoothing of B-spline Curves 

In order to fair or smooth at a particular spline knot, the method used 
must be local in nature. That is, when the method is applied, it only affects 
the curve in a small region surrounding the point which was smoothed. The 
technique used here is one proposed by Farin [6], [7]. First, the local fairness 
is calculated at each knot in the spline; then the knot with the largest value of 
e is chosen as the knot at which smoothing will take place. Let this knot be 
the knot associated with the B-spline control point dj. The knot is then re- 
moved from the knot sequence and a new location dj for the control point dj is 
calculated. Then the knot is reinserted into the knot sequence so that the 
number of spline segments remains the same as in the original curve. The 
criterion for the selection of a new location for the control point dj is third de- 
rivative continuity at the knot. Mathematically, this amounts to equating the 



left and right hand limits of the third derivative at the spline knot ( e.g. driv- 
ing e to zero ). The left and right hand limits of third derivative of a B-spline 
curve at the ith knot are given by 


p r u (°) = 77. 6 _ „ TtKi**/— i + "Hi d i + "H^i+1 + m AA d i + il 

\ U i + \ 11 i) 

mo) = yiKl^-2 + "*42*1-1 + "U^i + m AA d l + \\ ■ 

\ u i u i- V 

The values of coefficients of the d terms are calculated by equation (2.60). The 
new location of the control point di is given by equation (4.8): 

d = fr ? ill Z u ^ l ‘ + ~ U i-J r i (4.8) 

' (u i+2 -Ui_ 2 ) 

with the points li and rj given by 

I = (“/ + ! ~ ~ (“/ + 2 ~ "i-lfli-l 

' (", ~ «i- 3) 

= K+i ~ “/-3H+1 ~ (“,-+i ~ “/K +2 
(«i - u i - 3 ) 

Figure 4.1 illustrates the geometry behind equation (4.8). 

The curve smoothing routine developed for EAGLE takes, as input, the 
control vertices and the knot vector which define a B-spline curve. The rou- 
tine then smooths the spline as directed. The user may specify the number of 
smoothing passes to make, as well as a movement tolerance for points on the 
curve. The user may also specify whether or not the routine is to interrogate 
the curve and smooth only at points with large discontinuities. If the entire 
curve is to be smoothed, the user instructs the routine accordingly and all ver- 
tices on the spline, with the exception of the first two, and last two are 
smoothed according to equation (4.8). 
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Figure 4.1 Geometric interpretation of the B-spline smoother. 

The movement tolerance, mentioned earlier, enables the user to specify 
the maximum distance which any point on the curve is allowed to move. This 
option may be used where geometry perturbations must be kept within certain 
limits. Since the smoother moves the control vertices of the spline, maximum 
allowable movement of a point on the curve must be translated to an allowable 
perturbation limit for the curves control vertices. This is accomplished by use 
of equation (4.9): 



where 


(4.9) 
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Here, ej is the displacement vector for a control vertex necessary to displace a 
point on the curve by the vector 6xj. Using equation (4.9), a maximum dis- 
placement is calculated for each control vertex on the spline. The movement of 
the control vertices is monitored during the smoothing process to ensure that 
none of the vertices are displaced more than the allowable tolerance. If equa- 
tion (4.9) results in moving a control point farther than is allowed by the pre- 
scribed tolerance, the point is moved in the desired direction, but only as far as 
the tolerance allows. 

Smoothing of Tensor Product B-spline Surfaces 

Fairing of surfaces is accomplished by using a tensor product method. 
This amounts to using the curve smoothing routine to smooth the curve net 
that defines the surface. Each curve of constant v is smoothed and the results 
stored. Then, using the result of the first smoothing pass, each curve of con- 
stant u is smoothed. This procedure constitutes one smoothing pass. The re- 
sulting surface will be smoother than the original. The tensor product method 
is discussed in detail in Ref [7], and its application to surface smoothing is dis- 
cussed in Refs [6] and [7]. 

The surface smoothing routine developed for EAGLE takes, as input, 
the knot vectors and control vertices which define a tensor product B-spline 
surface. The surface is then smoothed using the tensor product method de- 
scribed above. The user may specify that only certain regions of the surface 
are to be smoothed. A region is defined by giving the indices of the vertices on 
two opposing comers of the patch. If the patch lies on the interior of the sur- 
face, the entire patch is smoothed. However, if one or more of the edges of the 
patch are coincident with one or more edges of the surface, the first two ver- 
tices from the edge of the surface are not smoothed. The user may, as with 


curves, specify a maximum allowable displacement. However, for surfaces, 
this value is the maximum allowable movement for a control vertex, not for a 
point which lies on the surface. 


CHAPTER V 


RESULTS AND CONCLUSIONS 

The material in this chapter is concerned with results obtained using 
the routines outlined in Chapter IV. Several examples of both curve and sur- 
face smoothing are provided as well as the results obtained from the respec- 
tive interrogation algorithms. Additional information is provided by results 
obtained from a numerical simulation effected on a NACA 0012 airfoil both 
before and after the smoothing algorithm had been applied. 

Curve Smoothing Results 
Waisted Body 

The waisted body represents a simple geometric entity described by 
four polynomial segments. Data for this geometry were given as a set of 120 
raw data points obtained from evaluation of the piecewise polynomial function 
describing the curve. The data points were then splined using a non-uniform 
B-spline; this resulted in a spline curve consisting of 120 control vertices and 
the associated knot vector. The curve smoothing routine was applied and the 
results are presented in Figure 6.1 through Figure 5.4. 

Figure 6.1 shows the control polygon of the spline both before and after 
smoothing. As evidenced by Figure 6.1, the geometry is changed very little by 
the smoothing process, however, Figure 6.2 illustrates the change in the 
splines curvature plot. Note how slope discontinuities in the curvature plot 
are virtually eliminated by the smoothing process. Figure 6.3 demonstrates 
the effect of smoothing on the metric coefficient r}y The degree by which the 
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geometry was perturbed is indicated by Figure 5.4. Note that the maximum 
point displacement was only about 1.6 percent of the total arclength of the 
curve. If Figure 6.1 is examined closely, it is hard to distinguish ( visually ) 
any change in the geometry due to smoothing; however, Figure 6.2 and Figure 
6.3 show that the continuity of the spline was substantially improved. 

Digitized Curve 

Figure 6.6 is a spline curve which was obtained from points digitized 
from an physical model Examination of the curve reveals visible discontinui- 
ties in first derivative. A better indication of the splines condition is given by 
the curvature plots in Figures 6.7 and 6.9. 

This data set was smoothed in two different ways. First, the smoother 
was allowed to pick the point with the largest third derivative discontinuity 
and smooth there. The results after ten passes using interrogation and point 
selection are shown in Figures 6.6. and 6.7. The actual spline curve is shown 
in Figure 6.6, while curvature plots before and after smoothing are shown in 
Figure 6.7. For the second case, the curve was smoothed without using the 
interrogation algorithm. In both cases the same number of smoothing passes 
were used. Figure 6.8 shows the spline after smoothing while Figure 6.9 
shows plots of curvature both before and after smoothing. 

NACA 0012 Airfoil 

The purpose of including these routines in a numerical grid generation 
program is to help enhance the quality of the geometry model used for numeri- 
cal simulation. Therefore, it is appropriate that some indication be given of 
how smoothing affects the results obtained from simulation programs. With 
this purpose in mind, data for a NACA 0012 airfoil were taken from Ref[24] 
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and splined. The data set for the airfoil consisted of 121 data points generated 
from the spline. The spline was then smoothed and another 121 data points 
generated. The smoothed and unsmoothed data sets are presented in Figure 
6.10. Again, close examination is required in order to. visually detect changes 
in the geometry ( see Figure 6.11 ). Figure 6.12 shows curvatures for the 
smoothed and unsmoothed geometries. 

To test how smoothing effects numerical simulations, an incompressible 
flow solution was run on both geometries. The development of the code used 
in to obtain the solution is documented in Ref[25]. Plots of pressure coefficient 
(CP) vs. X/C at zero degrees angle of attack for the two cases are shown in Fig- 
ures 6.13 and 6.14. Notice the smooth variation of CP for the modified geome- 
try. Also note that the CP curve for the smoothed airfoil still agrees with ex- 
perimental data. Color contour plots of the pressure field for the two cases are 
given in Figures 6.16 and 6.16. The simulation was also run at five degrees 
angle of attack. Plots of the pressure field near the leading edge of the wing 
section are shown in Figures 6.17 and 6.18. Notice that, for the unsmoothed 
geometry, the low pressure region on the upper surface exhibits ripples, 
whereas the same region on the smoothed geometry does not. 

Surface Smoothing Results 
Flat Plate 

A flat plate represents the simplest test for the surface smoother. Fig- 
ures 6.19 through 6.22 show a plane to which perturbations have been ap- 
plied. The surface in Figure 6.19 appears to be smooth and free of geometric 
irregularities. However, examination of its absolute curvature, shown in Fig- 
ure 6.20, reveals that the surface is actually rippled. After applying the fair- 
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ing algorithm, the new surface, shown in Figure 6.21 appears unchanged, but 
the curvature plot in Figure 5.22 shows that the irregularities have been re- 
moved. 


Sculpted Surface 

Figure 6.23 shows a surface patch digitized from a model of a lifting 
body vehicle. The surface is fairly smooth as evidenced by the plot of absolute 
curvature in Figure 6.24. The geometry was smoothed and the results shown 
in Figures 6.25 and 6.26. Additional options available for surface interroga- 
tion are shown by plots of first principal curvature and Gaussian curvature in 
Figures 6.27 and 6.28 respectively. 

F-16 Aft Quarter 

The last example uses a digitized surface from the aft section of an F-15 
fighter aircraft. The original surface is plotted in Figure 6.29 and contour 
plots of Gaussian and absolute curvature are shown in Figures 6.30 and 5.31 
respectively. The smoothed surface is shown in Figure 6.32 along with plots of 
Gaussian and absolute curvature in Figures 6.33 and 6.34 respectively. 

Conclusions and Recommendations 

The interrogation and smoothing routines, outlined herein, represent a 
powerful tool for evaluation of curve and surface quality as well as integrity 
enhancement of geometry data. These routines are limited, however, since 
they are unable to account for intended discontinuities. Research is still need- 
ed in order to produce an automated method for improving the quality of curve 
and surface representations without destroying desired geometric irregulari- 
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ties. At present, the only viable method is the use of these routines in an in- 
teractive environment. 
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Curvature of Smoothed Curve 
Curvature of input Curve 



Curvature plots for unsmoothed data and data smoothed using the interrogation algorithm 








Results of smoothing without using interrogation and point picking 





Curvature of Smoothed Data 
Curvature of Input Data 



Curvature plots for original data and data smoothed without interrogation and picking 





Before Smoothing 



B-spline control polygon for original and smoothed NACA 0012 airfoil 











Curvature of Smoothed Data 
Curvature of Input Data 



Curvature plots for the original data and smoothed data from the NACA 0012 wing section 







plot for original NACA 0012 data 







plot for the smoothed NACA 0012 data 
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ORIGINAL PAGE IS 
OF POOAQUAUTY 


Figure 5.15 Pressure field near the leading edge of the original NACA 0012 airfoil at 
zero degrees A.O.A. 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


Figure 5.17 Pressure field near the leading edge of the original NACA 0012 airfoil at five degrees A.O.A. 
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OF POOK QUALITY 


Figure 5.18 Pressure field necgr the leading edge of the smoothed NACA 0012 airfoil at five degrees A.O.A. 
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Figure 6.22 Color contours of absolute curvature on the smoothed planar surface 
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Figure 5.23 Original surface data for the lifting body surface patch 
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Figure 5.24 Absolute curvature on the original lifting body surface spline 
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Figure 6.29 Original surface data from the digitized F-15 aft quarter 
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Figure 6.31 Absolute curvature on the original F-16 aft quarter surface spline 
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Figure 5.34 Absolute curvature on the smoothed F-15 aft quarter surface spline 
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